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Abstract. Selfconsistent magnetospheres of rotating cosmic magnets (neutron 
stars/pulsars) with arbitrary inclination of the magnetic against the rotation axis are 
considered. Present studies concentrate on the regime dominated by the force-free 
surface (FFS). A macroscopic fluid description is applied and radiation reaction is 
taken into account. As in earlier work of our group, a 'standard set of parameters' is 
used. Under these conditions, the following features are found among other results: 
global charge separation exists for all degrees of inclination of the magnetic against the 
rotation axis; clouds of different charge are seperated by regions of vanishing particle 
number density; as expected, test particles inserted into the latter regions propagate 
into one of the adjacent clouds; strong polodial currents exist; locally averaged particle 
energies for protons typically range up to 10^^ — 10^^ eV, depending on the angle of 
inclination. 



PACS numbers: 97.60.Jd Neutron stars, 98.70.Sa Cosmic rays, 52.60.+h Relativistic 
plasma, 52.65.-y Plasma simulation 52.25.Wz Nonneutral Plasmas, 
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1. Introduction 

Since the discovery of pulsars in 1968 [1] and their interpretation as a rapidly rotating 
magnetized neutron stars in the same year [2, 3], these compact objects are under 
discussion as powerful accelerators of ultra high energy cosmic ray particles. 

Prom the very beginning of neutron star physics, work on the dynamics of 
electrically charged particles accelerated in the corresponding electromagnetic fields 
proceeded on two stages: (1) test particle dynamics in the vacuum fields of rotating 
magnets, e.g. [4, 5], demonstrating fundamental mechanisms and (2) self consistent 
plasma dynamics, e.g. [6], reproducing certain aspects of the structure and evolution of 
neutron star magnetosphcres. 

Since then, numerous papers have been published, investigating these matters in 
great detail, which we will not be able to discuss or even just to mention in this 
introduction to our present paper. Here we shall concentrate on stage (2), on relativistic 
plasma dynamics in a regime governed by the force-free surface of a homogeneously 
magnetized, rapidly rotating sphere with parameters typical for neutron stars. 

The notation of a force-free surface (FFS) refers to the dynamic of an electrically 
charged (test) particle within given electromagnetic fields. By definition the FFS is 
generated by those points of configuration space at which the Lorcntz-force acting on 
that particle through given electromagnetic fields vanishes. While in published literature 
a particle that happens to be at such a point often referred to as 'force-free' (ff), we prefer 
- in view of the presence of other types of electromagnetic forces (radiation reaction 
forces) - to speak of a Lorentz-force-free (Lff) particle in that situation. If B and 
E 7^ J, as is the case in fields considered here, a particle is Lff for E + [/3, B] = 0, 
i.e., for (I) (E,B) = and (II) (E,^) = and (III) ([E,H],/3) = [E^]. E is the 
electric field vector, B the vector of the magnetic induction and /3 is the velocity in 
units of the velocity of light. For |E| <^ |B|, as is the case under premises adopted 
here, one may expect the second and the third conditions to be inherently fulfilled to 
some approximation so that the FFS then is caracterized solely by the first condition, 
(E,B) = 0. 

From the early works of [7, 8] magnets rotating in the vacuum with the vector of 
magnetic dipole moment inclined against their respective rotation axis are known to 
create such FFS, of which some segments can act as particle traps and thus may have 
strong bearings on the formation of a neutron star magnetosphere, at least within a 
certain range of distance from its surface. 

In what follow, it will be useful to distinguish the special case of aligned rotators, 
i.e. rotating magnets with the magnetic axis parallel to the rotation axis {parallel 
rotators) or antiparallel to the rotation axis {antiparallel rotators), from the general 
case of inclined rotators and from rotating magnets with the magnetic axis orthogonal 
to the rotation axis {orthogonal rotators). On stage (2), a considerable number of 

I Boldface letters are for vectors in 3-dimensional euklidean space. (•,•) denotes the scalar product, 
[•, •] the vector product 
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investigations has been published on ahgned rotators. Some of these papers will be 
mentioned below. But only few are devoted to inclined or orthogonal rotators. In 
the latter case, obviously, considerable formal and numerical complications arise from 
the lack of rotational symmetry. Analytical approaches were used, for example, by 
[8, 9, 10, 11]. 

In the special case of aligned rotators, due to axial symmetry, theoretical 
results are achieved much easier, even on stage (2). Analytical methods have been 
applied to the structure of the magnetosphere of aligned rotators for example by 
[6, 7, 12, 13, 14, 15, 16, 17, 18]. Also, numerical studies on that matter have been 
performed by [19, 20, 21, 22, 23, 24]. 

On stage (1) in an earlier work of our group [25] the authors have integrated 
numerically the equation of motion for individual test particles within the regime of 
the FFS and with no restrictions on the relative orientation of the rotational towards 
the magnetic axis. Results confirm that velocity components orthogonal to the magnetic 
field vector are efficiently damped by radiation reaction. Thus, test particles tend to 
follow magnetic field lines, as suggested earlier by [7]. In a certain class of orbits they 
oscillate about the FFS, while moving along magnetic field lines. The amplitude of these 
oscillations decreases through radiation losses. Ultimately, in the subsequent regime of 
lower energy (and on a much larger time scale) particles become subject to drift in 
azimuthal direction. 

On stage (2) in a second paper of our group [26] the authors have introduced 
a numerical iterative approach to reproduce sequences of quasi stable plasma 
configurations forming under the influence of the FFS. In iterative steps charged particles 
were allowed to be ejected from surface elements of the rotating sphere in quantities 
locally proportional to the magnitude of the electric vector component normal to the 
respective surface element. These particles were then allowed to move freely along the 
appropriate magnetic field lines and to settle down where the projection of the electric 
onto the magnetic vector vanishes. Thus, without making use of the equation of motion, 
co-rotating, quasi stable, charge separated clouds were reproduced, in consistency with 
earlier results mentioned above. 

In our present work we proceed one step further taking into account in a 
selfconsistent way virtually all effects of relatwistic particle dynamics, including 
radiation reaction and effects of special relativity as, for example, retardation. Our 
numerial approach is designed to describe the evolution of locally averaged particle 
densities, since velocity dispersion is not taken into account. Again, a magnetosphere is 
allowed to build up from the initial vacuum through particle ejection from the spherical 
surface, similar to the procedure described above. Here particles are allowed to be 
ejected from surface elements, given an appropriate direction of the electric vector, in 
quantities locally proportional to the magnitude of the electric vector component 
projected onto the magnetic field line. 

Thereby, we intend to clearify, on stage (2), the evolution, selfconsistent structure 
and stabihty properties of plasma configurations forming within the regime of the FFS 
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of a rotating cosmic magnet with no restrictions on the relative orientation of the 
rotational towards the magnetic axis. Also, we want to evaluate mean energy values 
locally achieved by particles in that regime. 

As in earlier works of our group, e.g. [27, 28, 29], we apply a model represented 
by a rotating, ideally conducting, homogeneously magnetized sphere with arbitrary 
inclination of the magnetic against the rotation axis. A 'standard set of parameters' 
representing well-known properties of typical neutron stars is attributed to this model: 
the stellar mass which is taken equal to the solar mass rriN — rusmn the stellar radius 
rjv = lO^cm, the angular velocity uj — 207rs~"^, and the magnetic dipole moment 
^ = 1030 G cm^ §. 

In order to investigate the regime of the FFS with an appropriate resolution we 
concentrate on the near zone of the neutron star up to 20r7v. 

From preceeding estimates as well as from subsequent simulations gravitational 
forces exerted by the rotating neutron star and by the magnetosphere itself onto 
individual particles, as well as effects of general relativity were found to be negligible 
for the standard set of parameters. Contributions to the Lorentz-force originating from 
magnetic fields created by magnetospheric particle currents can also be neglected, in 
agreement with earlier conclusions of [22]. Spontaneous pair creation still turns out to 
be insignificant, even within the very strong electromagnetic fields of polar regions. 

In chapter 2 of what follows, we rediscuss force-free surfaces associated with 
rotating magnets. Thereafter, in chapter 3, we display equations of individual as well as 
of collective particle motion in terms of a non-neutral two-component fluid description, 
and we then proceed to a description of appropriate tools for numerical treatment in 
chapter 4. Results will be given in chapter 5 and subsequently discussed in chapter 6. 

2. Vacuum Fields and Force-Free Surfaces 

The vacuum solution of Maxwell's equations for a homogeneously magnetized (ideally 
conducting) sphere, rotating with its vector of angular velocity uj inclined relative to its 
vector of magnetic dipol moment ^ by the angle as evaluated in [30], called Deutsch- 
field, may be applied here in the near-field approximation. In addition, to account for 
the global electric charge of the rotating sphere, an electric monopol contribution is 
introduced. For Details about the Deutsch-field in the near-field approximation with a 
global electric charge of the rotating sphere we refer to [26]. 

In the case of an ideally conducting sphere, the topography of the exterior field 
is known to be independent of the form of interior magnetization. An interior central 
magnetic point dipole, which may be chosen as an alternative model of magnetization, 

§ Gaussian units are used throughout this paper. Thus, electric and magnetic field strengths are 
measured in units of IG = 300V/cm. For the standard set of parameters, the magnetic field strength 
is about Bp ~ 2 • 10^^ G and the electricic field strength is approximately Ep ~ 10^°G in the polar 
region. Under given parameter values, the radius of the light-cylinder (often referred to as 'light 
radius') n, = w /c, outside of which corotation cannot exist, is ri, = 5000 km, corresponding to almost 
the radius of the earth, for comparison. 
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Figure 1. Force-Free Surface. Left: X = 60°, Qs = h/vl- Middle: X = 60°, 
qs = 0.3 iJ./rL- Right: X — 90°, = /x/rL. The unit of length on both axes: rjv- 



obviously would create the same electromagnetic vacuum field configuration. But the 
electric surface charge as well as the interior charge evoked by rotation clearly depend 
on the form of magnetization. 

A central interior magnetic point dipole, for example, is consistent with the 
total interior electric charge Qi — |^cosx and with the surface charge density 
[cos X cos^ ?9 + sinxcos^sini?cos(</' — corresponding to the total 



surface charge Qs = — J^cosx ||. 

Alternatively, for a homogeneously magnetized sphere, as considered throughout 



this paper, the surface charge density is a 
5 sin X cos sin i) cos(v^ — to t)] and 



;3-[cOSx(5 COS^ — 3) + 



iil COS X = -Qs ■ 

The discontinuity of the tangential component of the magnetic surface field creates 
an electric surface current which results negligible under conditions given here. 

With the projection of the electric onto the magnetic vector, written in the /i- 
system, 

Hk^ 1 



(E,B) 
B 



the FFS is given by 

r 



(u ^4 3-T '<^I^n) (sin X sin cos a -2 cos V') 

[kr)^ vl + ocos"' ip L 

+ 4(cos X cos — sin X sin cos A) cos^ ip — sin x sin cos A] 



(sin X sin ip cos A — 4(cos x cos ip sin x cos A sin ip) cos^ ip) 
(sin X sin 'ip cos A — cos ip) 



where := qs'^ is the dimensionless form of the total electric charge of the sphere. 



Some examples for the chape of the FSS are illustrated in figure 1. 

II Here we make use of two sets of spherical coordinates: one is referred to as the w-system (r, 'd, 'P), 
where r is the radial coordinate, is the angle measured against the rotation axis, and is the angle 
relative to the plane spanned by the Xq and yo axis, at rest in a chosen inertial frame of reference. 

The other set of spherical coordinates is referred to as the /x-system (r, tp,X), in which ^ is the angle 
relative to the magnetic dipol axis, and A is the angle against the plane spanned by the fi and u> axis. 
As a consequence of these definitions, A = ^ — u}[t— (r — rjv)/c]. 
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3. Non— Neutral Two— Component Plasma Fluid 

3.1. Equations of Motion for the Plasma 

We consider a two-component ideal fluid, whereby each component (identified by the 
index s = 1, 2)consists of N identical (classical) particles, electrons and protons in 
this case. For a relativistic macroscopic description of fluid motion in terms of the 
one-particle distribution function Ts{x^-iP^) we start with the one-particle Liouville- 
equation 

where = ct is the time coordinate, are the spatial coordinates, p'* are the four 
corresponding components of momentum Using the covariant form of Hamilton's 
equations, dx^/dr — dH/dpn and dp^/dr — —dH/dx^, where H is the Hamilton- 
function, one may write the covariant Vlasov-equation in the form 

^ (•F.x^) + #-(^.p'') = 0. (3) 



With the definition of the one-particle distribution function fs{x'^,p^) '■ = 
f'^oo ^s{x^-,P^) dp^ the Vlasov-equation results (with the assumption ^ ^ for 



p'^ oo) in: 

If = dx^/dr are the four components of velocity, = du^/dr those of the 
acceleration, while 7 is the Lorentz-factor and r the proper time, the Vlasov-equation 
can be represented by 

|(T/,)+cA(«-/.) + A(„./.) = o. (5) 

The distribution function fs{x^,u^) delivers by integration the average number density 

ns{x^) = j fs{x^,u')d^u (6) 

and from there the average, macroscopic (four-component) velocity vector and the 
electric charge density is given by: 

<(x'^) = n;' j u'^ /,(x^ u') d'u , ^(x^) = 5^ risix'^). (7) 

s 

Likewise, wc make use of the average, macroscopic (three-component) velocity vector 
vl{x'^) to define the electric current density 

f{xn ^^ej' v^fsd^u = ^sv\{x^) (8) 

s ^ s 

^ The signature of the metric tensor is (1,-1,-1,-1). Throughout this paper, Greek indices are running 

from to 3, latin indices from 1 to 3. As stated before; efFc;cts of general relativity (including gravitation) 
are left unaccounted for since, with the standard set of parameters the Schwarzschild-radius is only 
about Ts ~ 0.3 rjv- 
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constituting thereby the four components of the electric current density vector field 

jf^ = {cQ,f). 

The latter, in consistency with (5), act as source terms in Maxwell's-equations 
Q^pc,t3 ^ ^Q^p*c,fi ^ Q ^ ^Yieie F"^ = d'^A'^ - d'^A'' is the electromagnetic field 
tensor and F*°'^ = ^e"'^^^ F^s its dual counterpart. A^ represents the electromagnetic 
potential and e°'^'^^ the Levi-Civita-tensor. 

The average number density ns{x^) and velocity ul{x^) fields of each of the two 
constituents of the plasma serve to describe macroscopic fiuid motion with the help of 
the continuity equation 

^""^ + ^(nsvl) = - / #-(a7.) d'u = Q, (9) 



dt dx^ J du 

obtained from (5) through integration over velocity space, adopting |a|/s for 
|u| oo. The energy-momentum equation 

^ + vi^n', = {-d,p^ + IC) , (10) 

ot ox^ [msns) ^ ' 

is obtained from (5) through multiplication with and integration over the velocity 
space, where 

p''i{x'^) = j ms{u' - ul){u^ - ui)fsd^u^msj u'u^ fgd^u - msUsului (11) 

is the pressure tensor and the rest mass of particles constituting the respective 
plasma component. 

In what follows, we assume a coUisionless, dispersionless (i.e. cold), relativistic 
plasma. This assumption appears not implausible, since with a neutron star surface 
temperature of about 10^ — 10^ K and a particle density of about ^ 10^^ cm~^, 
as suggested by the Goldreich & Julian-model (1969), the plasma parameter A = 
(47r/3)nsAD (where Az) = ^JUbT / (47rnse^) is the Debye-length) results in A a; 10'', 
so that A 1 is given. 

Thus, the right side of equation (10) reduces to the (volume) force term 

r^msja'fsd%, (12) 

to which the Lorentz-forcc is expected to deliver a major contribution JC^Lorentz — 
TjoF^'^u,^, where 770 = e/{rnc). 

3.2. Equations of Motion for Individual Electrically Charged Particles 

In addition to the Lorentz-force the classical equation of motion for an electrically 

charged relativistic particle subject to given ('external') electromagnetic fields has to 
account for radiation reaction forces. One such equation which frequently is called 
'Lorentz-Dirac (LD) equation' [31] (though Abraham-Lorentz (AL) equation would be 
historically more correct) may be given the form 
dii^ 

— = r]oF^^''u, + roG''''u, (13) 
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suggested by one of us [32] . Here, 



1^ 



(14) 



dr2 ' 

is the radiation force tensor and Tq = 2e^/ (3mc^) is the radiation constant. 
Unfortunately, this equation of motion exhibits serious deficiencies which have 
extensively been discussed in published literature and can be avoided through 
replacement of (14) by its first iteration often referred to as 'Lorentz-Dirac-Landau 
(LDL) equation' (due to its extensive discussion in [33] well known textbook). In this 
approximation the radiation tensor may be written 



(15) 



with -u^^ 



In our numerical work presented further below we adopt 



that locally (i.e. within appropriately small intervals of space and time coordinates) the 
electromagnetic field is (approximately) homogenous in space and constant in time, in 
which case the radiation tensor (15) further reduces to 

(16) 



— '^const — ^oVo 



F%F\u- 



+ F\F/u^u^u'' 



Under these premises individual particles of each of the two plasma components, 
according to s = 1 or 2, inside the proper interval of space and time coordinates 
(i.e. inside the corresponding 'volume element') are subject to the same electromagnetic 
forces. Consequently, the (macroscopic field) equations for the plasma fluid * 
immediately follows from (10) and (13) with (16) 



up, 



sd^u': = r}oF''%,s + Tovl 



+ F\F/uiu^,,u^ 



(17) 



3.3. Exact Solutions of the Equation of Motion for Individual Particles in Homogenous 
and Constant Fields 

In what follows we make use of exact analytical solutions of the equation of motion 
(16) for individual electrically charged particles in locally constant and homogenous 
electromagnetic fields [29]. Given a coordinate system with 63 = B, 62 = [E, B] and 
61 = [e2, 63] and excluding null fields jj (i.e. fields with simultaneously vanishing Lorentz- 
invariants (E, B) and — B^) and restricting further to 7^ and -B 7^ 0, the solution 
of (16) is given by: 

/ cosh Xt\ ( 



7 a(r) < 



\ 



13 j^sinh A t 
— P cosh A T 
-sinh A T 



+ C3 



V 



sinh A t\ 
P j^cosh A T 
— (5 sinh A t 
-cosh A T 



One of us [29] has argued that in a quantum-mechanical frame self-consistency of classical 
electrodynamics suggests (15) to be the correct form of the radiadion tensor G^" 

* Density effects on radiation and radiation reaction remain unaccounted for. Remarkably, since in (17) 
the pressure term is absent, the number density does not appear. Clearly, (17) governs the macroscopic 

velocity fields u^. 

ji Inside the magnetosphere, |E| <C |B| 
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+ 7 6(r) < 



f3 sin UT^ 

COS LOT 

— sin UJT 



+ C2 



f2 



P COS cjr\ 
sina;T 



cosu;r 



(18) 



E2+b2_^(52_e2)2+4 (E,B)2 

2EiB • 



2^ 



With A = sign(^3)£:\/-K^'-^') + W(^'-^T + 4(E,B)2, ^ ^ ^ 



e 

mc 



E,B) 



C2 



+ 1 , 
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V 



/?«0(0) + «2(o) 

/ 



and the radiation parameters a(r) = [{Cq — C|) — {Cf + C|) exp[— 2ro(A^ +uj'^)t]]~^ , 
6(r) = a(T) exp[-To(A^ +0;^) r], where C^tC"" = 1 . 

If the electric field is parallel (or antiparallel)to the magnetic one {Ei — 0) this 
solution reduces to: 

/ coshAT\ 



V sinh Xt J 

I \ 
sina;T 
cosa;T 



u^{t) = a(r) < 



u\0) 



I ) 



+ b{T) < 



AO) 



V 



+ u%0) 



I 



/ M 

cosa;T 
— sincur 

V 

/ sinhAT\ ' 



V cosh A rj ^ 



(19) 



The dynamics of a charged particle starting at rest is characterized by the following 
properties: Within a regime of small proper time, r <^ l/ro(cj^ + A^), the radiation 
parameters are about air) ^ b{r) ^ 1, while particle motion is a composit of gyration 
around magnetic field lines, [E, B] -drift, and acceleration along magnetic field lines, due 
to the projection of the electric field vector onto to the tangent to the magnetic field 
line. 

Within a regime of large proper time, r — > 00, gyrations around magnetic field 
lines are damped corresponding to 6(r) — > 0, while acceleration along magnetic field 
lines proceeds correspondig to a(T) 1, superimposed by [E,B]-drift and curvature 
drift. 

In the magnetospherical regime considered here (|E| -C |B|), where uu ~ ^l,X ~ 



f^L-i sin a and 



1 f 5 ■ 10 s for electrons 

ro(a;2 + A^) \ 3 • 10~^s for protons , 



(20) 
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gyrations around magnetic field lines are damped. 

4. Numerical Procedures on the Grid 

4.1. Velocity Components 

Evolution in time of the four components of the velocity vector u"{t) of individual 
electrically charged particles, from 'u"(to) to u°'{ti) with At = ti — (with substitution 
of proper time r by the time coordinate t, involving the integration of vP{t) over r) 
according to what was discussed before, is governed by the equation of motion (13) with 
(15) (within appropriately small intervals of space and time coordinates). The same 
holds for a plasma velocity field, in which, according to the dispersion-free fluid model 
used here, individual particles represent co-moving fluid elements, 

«"(Ti)-«"(ro)=M"(ti)-M"(to)= / —^dt= (dt + v,it,x)&)u''it,ic)dt.{21) 

Jto at Jto ^ ' 

The first term in the last integral of (21) refers to the explicit time dependence of the 

velocity field (i.e. the change with time at a fixed point x in space), whereas the second 

term describes its implicit time dependence (i.e. the additional change with time of a 

co-mpoving particular fiuid element). 

Here we are interested in the change in time of the velocity field at a fixed point 

X in configuration space. Due to the given initial condition ^"(ro = 0) = M"(to,x) in 

each iteration step (where u"(t,x) denotes the velocity at a given point in configuration 

space) it follows: 

M"(tl, X) = M"(ti) - Vj{t, x)^' M"(t, x) dt . (22) 
Jto 

Velocity components are required on grid points of a spherical grid, where the 
components of the electric field vector are suggested to be known from electric charge 
density averaged over each individual cell. 

The integral in (22) can be evaluated in a first order scheme in time with 
V := -{ti-t) dtv = 1, 

ti°(ii,x) = u^{t,) - {ti-to)vj{to,^)d^u''{to,^) - f\ti-t) dt[v^{t,^)&u''{t,^)] dt (23) 

Jto 

by neglecting the integral in this evaluation. Differentiation with respect to spatial 
coordinates is performed numerically in a second order centered difference scheme 
(boundaries are treated separately, also in second order) . The numerical scheme results 
in 

ii"(ii,x,) = u^itr)- (24) 

At / 1 1 1 \ 1 . ^ , .a, 

7 V Ar rAi9 rsmi^A'/'/ 2 

where the index i now characterizes the i*'* grid point (e.g. with respect to Cr-direction, 
etc.). 
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4-2. Electromagnetic Field Components 

Numerical integration of electromagnetic field equations is performed applying a scheme 
developped in our group by [34], implying a (complete and orthonormal) system of 
spherical vector harmonics Pnm,Bnin, Cnm (e.g. [35]), 

^n(n + l) 

CnM, <P)^ - [Br, Br,m{^, <p)] = , / [V, rX^J , 

^Jn[n + 1) 

where n e {0, 1, 2, ...}, m e {—n, 0, n} and the spherical harmonics X^mi'^, f) ^-re 
defined with the help of the associated Legendre-function Xnm{,'d-iV) — e*"**^ P^(cosi9) 
permitting the expansion of the elctromagnetic vector potential 

A(r, 1?, ^,t)^Y. [Pnm{r, i)Pnm(l?, ^) + &nm(r, i)B„^(^?, + c„^(r, t)C„^(i?, V') ] . (25) 



For example, the vector potential of the Deutsch-field can be represented in the form 
A(r, ^M-t'^ cos X P20 - ^e-^* sin x P21 



+ COSXB20 -^e '*^-FF7 — rsmxBai 

+ lcosxCio + e-*^^i^sinxCn, (26) 

with /ii(r) = -e-^ , h2{r) = ^ e-^i±2^^$^ , //^(r) = e'r {^r-r^)+f-^r-) 

Prom there the magnetic field vector is calculated from B = [V, A]. Exploiting 

the gauge invariance of the four component vector potential to eleminate the electric 

field vector is calculated by E = — 9tA. 

To evaluate the total electromagnetic field in terms of its expansion coefficients, 

we add the expansion coefficients of the Deutsch (vacuum) field to the expansion 

coefficients of the the plasma (different gauges are used for the two components). 

The electric (scalar) potential AQ{r^il)^a) is determined from the charge density 

0{r,i3,'^)- Ao{r) = / jf^^'^'^^' > assuming Ao(r7v) = Ao{oo). Furthermore, 
= '}ln,mWnm-p^ Xl^X'^j' Xnrni'd,'^), wherc u)„„, = {u - m)\ / {u + m)\ and, by 

definition, r< (r>) refers to the lower (upper) limit of the considered range of |r| and 

|r'|, respectively, so that 

Ao{v) = Y.^nm{^,V) j Wnm^Q{v')Xl^{§',V')d'r' = Y.^nm{r)X^m{^,V) (27) 

n,m n,m 

with Xnm{r) = JWnm:^Q{r')X*^{'d',V')d^r' = J Qnmir') r'"^ dr' . Here ^„m(r') = 
I J WnmQ{'^')X*^{i)' ,^') sm{i^')di!}'dip' are the expansion coefficients of charge density. 
For a thin spherical shell of thickness Ar these expansion coefficients are 



Xnm{r) ^ Ar r| Qnm{rj) = Qnmirj) (28) 
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with qnm{rj) = Arr]gnm{rj)- 

The modes of image charges can be seen as additional modes of the surface charge 
density. The image charge at position ri is characterized by Vgp = rj^/ri and Q{rsp) = 
-{rN/ri)Q{ri). In general, r^p = r< so that qnm{rsp) = -^qnm{ri), Xnm{rsp) = 

3#Tq„„(r7v), with q„^(r7v) = --^qnm{ri)- 
^> 'i 

The total potential at grid points r = {ri^'dj,^^) inside the spherical volume is 
written Anmij'-i) = J2j<^irj'qnm{'>^j) and at grid points outside the spherical volume 
Bnm{ri) = J2j>i -ikrqnmirj). Prom that, the potential is 

3 

M^i: Vj: '&k) = '^nmin) X^mi^j, Vk) , 

n,m 

with T>nm = Anm'^i^"'^^'' + I3nm{fi)i"i', resulting in the coefficients of the electric field 

Pmn{r) = -drVnmir), hmn{r) = -^^^^P^I'nm(r), C^„(r) = 0. 

4.3. Continuity Equation 

In order to integrate the continuity equation dixj'^ = the method of flux corrected 
transport (FCT) is used. This method bases on an algorithm developped by [36, 37, 38]. 
In a ffist step a numerical 'low-order' scheme is applied introducing sufficient local 
numerical diffusion in order the get the numerical integration of the transport equation 
stable and monotonous. In a second 'high-order' step the introduced numerical diffusion 
is ehminated as far as possible. For details regarding FCT we refer to [38] . 

We extended this procedure to three dimensions, which is shown in Appendix A. At 
the end we get a stable conservative discretisation scheme with low numerical diffusion 
to integrate the continuity equation numerically. 

4.4- Particle Injection 

Of the three frequently discussed injection mechanisms - (1) emission from the spherical 
surface (as determined by the electric field topography at the surface and surface charge 
density), (2) invasion of particles from outer regions, and (3) electron-positron pair 
creation from photon decay - we need to consider only (1) here. The rate of particle 
injection from the surface is chosen to be proportional to the magnitude of the electric 
field component projected onto the tangent to the corresponding magnetic field line 
E\\ = sign(cos-?/')(E, B)/B, if the sign of surface charge density agrees with the sign of 
£■11 at the respective point on the surface. 

4-5. Reproducing the Magneto spheric Configuration 

All simulations carried out in order to get the results presented in this paper are started 
from the vacuum case. For X = (x = tt) the electric field in vacuum allows emission 
of electrons (protons) only, while emission regions of electrons and protons are equally 
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large for x = |. For < x < | the emission of electrons and for | < x < tt the emission 
of protons predominates. 

We study the magnetosphere of an initially non-charged, homogenous magnetized 
sphere up to 20r7v with the standard set of parameters. The following grid sizes are 
used in the numerical simulation: 

At' = 5.0 • 10"^ if X = 0, X = TT and At' = 2.5 • 10"^ if < x < tt 
Ar' = 3.99 • 10"^ 

A^?' = AV'' = 2.06 • 10"^ i6 and A^?' = A^^' = 4.12 • 10"^ ob, 

where ib destigates the inner border and ob the outer border. This resolution implies 
that e.g. for the (anti)parallel rotator a fluid element can cross the radial simulation 
extension 50 times during simulation time, which corresponds to a 114.6° rotation 
of the rotating oject. For the oblique rotator, the simulation time corresponds to a 
rotation of 57.3°, so that a fluid element can cross the radial simulation extension 25 
times. Consequently, the simulation time is large enough to reproduce existing quasi- 
stationary magnetospheres. We study the rotator exemplary for several inclination 
angles: x = 0°, 30°, 45°, 60°, 75°, 90°, 120°. 

5. Results 

We present and discuss our results under the aspect whether the formation of quasi- 
stationary magnetospheres can be verified and if so, what can be said of its structural 
features. In this context we investigate the verification of the predictions of the Goldreich 
& Julian-model and other authors in the special case of the aligned rotator. In the case of 
the inclined and orthogonal rotators (in which little is known from published literature) 
we investigate the structure of the quasi-stationary magnetospheres as well and stress 
in general the question concerning the typical particle number densities, currents and 
average particle energies inside the magnetospheres in the regime of the FFS. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.5 1 1.5 

t [a> ^] Inclination angle % 

Figure 2. Left fig.: Sphere charge depending on the time for different inclination 
angle. Right fig.: Sphere charge depending on the inclination angle in the case of 
quasi-stationary magnetospheres. 
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Figure 2 (left) demonstrates the development of the total electric charge of the 
sphere (over time) for various values of the inclination angle X) ranging from x = 
to X = ^/2- Obviously, an asymptotic value is reached in each of these cases within 
time intervals well below t ^ 0.2 o;"^, corresponding to a rotation of ^ 11.2°. This 
behaviour can be seen as a clear indication for the rapid formation of quasi-stationary 
magnetospheric configurations for all inclination angles. 

Furthermore, a prametrization of the electric monopole of the rotating sphere after 
building quasi-stationary magnetospheres in term of the inclination angle is possible 
and leads to (figure 2, right): Qs — |^cosx. 

Likewise, diagramms on the left of figure 3 show the development of the electric 
vector component parallel to the tangent to the respective magnetic field line over time 
on the stellar surface. Different curves, plotted against the polar angle, correspond to 
difi^ercnt values of the time coordinate. The three diagramms from the top to the bottom 
of figure 3 are for different values of the inclination angle ranging from X = 0° , % = 60° 
to X = 90°. In all these cases the electric vector component parallel to the tangent to the 
respective magnetic field line vanishes for all inclination angles within a time interval 
of about t 0.20;"^, which indicates the formation of quasi stable magnetospheres. 
Analogous conclusions can be drawn from diagramms on the right of figure 3 for the 
electric charge density on the surface of the sphere. 

For an analysis of spatial electric charge density inside the quasi-stationary 
magnetospheric particle distributions the corresponding even modes (normalized to r^) 
in terms of coefficients of charge density (sec chapter 4.2) as a function of the radial 
coordinate are shown in figure 4. Different curves are for different modes, while different 
diagramms indicate different inclination angles. In the special case of the aligned rotator 
the space charge density is described very well only by the quadropol mode n — 2,m — 
in agreement with the predictions of the Goldreich & Julian-model as well as of the more 
recent work by [24] ff . For increasing inclination angle x other modes gain more and 
more importance, especially those with m 7^ which are responsible for non-axially 
symmetric constibutions, as shown in figure 4. 

For a better vizualization of spatial electric charge distribution inside the quasi- 
stationary configuration, in figure 5 the number densities of the electron fluid (on the 
left) and of the proton fluid (on the right) within the plane spanned by the magnetic 
and the rotation axis are shown for different inclination angles. 

The structures of the resulting quasi-stationary magnetospheres are dominated 
by the force-free surfaces for all inclination angle and are completly charge seperated 
devided by regions of vanishing particle number density, often refer as 'vacuum gaps'. 

In the case of the inclined and orthogonal rotator, < x < f , electrons are collected 

ft In the Goldreich & Juhan model [6] a force free magnetosphere ((E,B) = 0) and a co-rotating 
plasma inside the light-cylinder are assumed. With E = —[P^°'^°,'B] and /3'*°™ = (r/r/,) sin i?e<^, the 
electric field outside the sphere results in: E = — ■^|p-(— sin^ er -|- 2 cos 1? sin e^). This field is caused 
by the charge density ^gj = —:;^^P2{cos'd) and an additional sphere-monopol by = The 
electric potential is given by ^0 = ^0"°"° + Ag"""^™ = | Jii( 1 - P2(cosi?) ). 
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Figure 3. Electric parallel field E'^^ := sign(cos V')(E, B)/S on the surface of the 

sphere (left plots) and the surface charge density (right plots) depending on the polar- 
angle relative to the rotation axis (A = 0). From the top toword the bottom:x = 0°, 
X = 60°, X = 90°. 



between the rotation and the magnetic axis, nearby and in the plane spanned by these 
axes. The fluid includes the rotation axis, except for |x — || 0. Protons are collected 
between the equator plane relative to the rotation axis and the equator plane relative 
to the magnetic axis, once again nearby and in the plane spanned by these axes. Given 
I < X < TT the sign of the particles changes. 

We found corotation (relative to the surface of the sphere) for all inclination 
angle. Particle number densities inside these clouds, for the standard set of parameters, 
typically range up to 10^^ cm^"^. 

For a discussion on currents and particle acceleration within these clouds normalized 
velocity fields of electrons and protons are considered in figure 6 for various inclination 
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Figure 4. Coeffients of the space charge density in the quasi-stationary case for 
different inclination angles. Upper rows: m = 0; lower rows: m ^ 0. 



angles. Furthermore the projection of the electric vector onto the tangent to the 
magnetic field lines, (E, B)/|B|, is shown. The direction of the velocity field correlates 
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X [t:J<r| X [r_N] 



Figure 5. Particle density [cm^^] of the electron fluid (left fig. ) and proton fluid (right 
fig. ) mapped over logio(£' '^"^^ + 1) after a rotation ~ 57.3° for different inclination 
angles. 

monotonously with the sign of these projection of the electric vector. The averaged 
Lorentz-f actors for electrons and protons are shown in figure 7. Average values are 
calculated for spherical shells (the radius of a given shell is mapped on the abscissa). 
Different curves indicates different values of the inclination angle. 
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In the case of the parallel (antiparallel) rotator, where electrons (protons) are 
collected around the rotation axis and protons (electrons) around the equator plane, 
a polodial (outward directed) current along the magnetic field lines exists, consisting of 
electrons (protons). In the equator plane we found protons (electrons) diffuse out of the 
simulation volume. In the context of the polodial current in the case of the antiparallel 
rotator averaged protons energies up to 10^^ eV are found. 

For inclined and orthogonal rotators (0 < x < |) outward directed currents 
consisting of electrons along the magnetic field lines are found. The angle range 
according the polar angle of these currents decreases with decreasing values of |x — f |- 
With inclination angles 30° < x ^ 90° closed currents consisting of protons are observed 
in the range of a few r^i starting and ending at the surface of the sphere. Given 
^ < X < vr the sign of the particles changes. In the case of the 120°-rotator, investigated 
as an example for inclined rotators with X > f and with currents consiting of protons, 
averaged proton energies up to 10^'' eV has been proven. 

The influence of the [E, B] -drift on the structure of the magnetosphere is increasing 
with decreasing |x — f |. Particles of both sign (electrons and protons) are streaming 
due to this force back to the surface of the sphere. 

6. Discussion 

The resulting quasi-stationary magnetospheres are not global force-free. In regions 
with vanishing particle densities the vacuum electromagnetic field is approximately 
undisturbed by the non-neutral plasma. In regions with high particle number densities 
the projection of the electric vector onto the tangent to the magnetic field vanishes 
nearly, but not complete. Small deviations from a force-free situation leads due to the 
extrem strong electromagnetic fields immediately to high particle energies. All in all, 
independent of the inclination angle, highly rclativistic plasmas are found, which leads 
to the necessity to take the radiation reaction into account. 

In the case of the parallel and antiparallel rotator, due to the fact that we did 
not found particles with low odcr moderate particle energies, the usage of the ultra- 
relativistic approximation of the Lorentz-Dirac-Landau-equation by [22] and [23] in 
their numerical calcualtion is justifiable, shown by our studies. Puthermore, we can 
confirm the existence of clouds with different charges, seperated by regions of vanishing 
particle number density (vacuum-gaps), also found by [21], [22], [23], [26] and [24]. 
All these works including ours found no global force- free magnetosphere. Force- free 
magnetospheres are often used as an assumption in analytical works. The structure 
of the quasi- stationary magnetospheres predicted by [7] with analytical models can be 
confirmed with our studies, especially the importance of the FFS. Nevertheless, the 
closed polodial currents proposed by [7] can not be proven by our studies due to the 
limited simulation volumen of 20 rjv- The density distribution of the plasma is described 
very well by only one quadrupole mode, confirming with [6] in their amplitude and power 
law. 
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Figure 6. Fig. left: Normalized velocity field of the electron fluid. Fig. middle: 
Normalized velocity field of the proton fluid. Fig. right: Electric parallel field 
E^^ = (E, B)/i3. li^iil > 10 is mapped to ±10 in order to show the change in the 
sign. Prom top to bottom: x = 0°, x = 60°, x = 90° 

In the case of the inchned and orthogonal rotator only few literature is published. 
The by [9] proposed approach solving the stationary Valsov-equation and the Maxwell- 
equations selfconsistently led to an inside the light- cylinder corotating, charge seperated 
magnet osphere. The charge separation independent of the inclination angle was also 
found by [26], which can be confirmed by our present studies. Basicly, the structure of 
the quasi-stationary magnetospheres in our calculations confirms with those found by 
[8] analytically. Furthermore, we can attest the dependency of the electric monopole of 
the rotating sphere of the declination angle, first proposed by [8]. Beside this, for the 
first time with our work statements on particle engeries, currents and drifts based on 
numerical, full dynamical studies are possible. Due to the limited simulation volume 
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Figure 7. Averaged Lorentz-factors per radial sphere depending on the radius for 
different inclination angles. 



we can not prove if outside the light-cylinder closed currents are formed (starting and 
ending on the surface of the sphere). The knowledge of the global current system is 
important due to the fact that the [E, B]-drift is important at long time scales. As a 
consequence of these force charged particles move back to the surface and may change 
the electric monopol of the sphere. The structure of the FFS and resulting from this 
the structure of the magenetosphere is depending on the electric monopol of the sphere. 

The high particle energies in the existing polodial currents arc relevant for induced 
pair production. In the present paper we were able to ignore the resulting e^-plasma. 
Nevertheless, relativistic currents in these plasmas may cause microscopic instabilities 
which could explain the non-thermal radiation in real pulsar magnetospheres. 
Analytical studies on this topic (e.g. [39]) assume Lorentz-factors in these currents 
of 7 10^ and particle densities in the scale of the Goldreich & Julian-density, as were 
confirmed by our studies. 

Regarding the discussion about neutron stars as cosmic accelerators for ultra high 
energy cosmic ray particles the very high particle energies proven by our studies are 
remarkably. I.e., in the case of the 120°-rotator we found averaged proton energies 
up to lO"*^^ cV. In the case of higher magnetic fields, than given by the 'standard set 
of parameters' higher particle energies are possible. But, investigating these cases one 
schould prove whether a classical approach is applicable. It is important to note, that 
due to the limited simulation volume used in our work we can not predict if these 
high energy particle are able to leave the neutron star magnetoshphere and if, at which 
particle energies. 

7. Summary 

In this paper we studied relativistic magnetospheres of rotating cosmic magnets (neutron 
stars/pulsars) with arbitrary inclination of the magnetic against the rotation axis. 
Concentrating on the regime dominated by the force-free surface (FFS) we developed 
a macroscopic description of a cold, coUisionless two-component fluid, consisting of 
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electrons and protons, taken into account the radiation reaction and carried out 
selfconsistent numerical calculation of relativistic magnetospheres of neutron stars. 

According the first two moments of the relativistic Vlasov-equation the equation of 
motion of the fiuid components are derived. Under the assumption of a cold, collisionless 
plasma considering the radiation reaction is possible. Due to missing velocity dispersion 
of the fluid components the radiation reaction term of the equation of motion for a single, 
charged particle can be added to the equation of motion of the fluids in the macroscopic 
description. Dealing with near zone of rotating cosmic objects up 20 sphere radii the 
influence of existing currents to the vaccum magnetic dipol field can be neglected. 

Beside the investigations regarding the parallel and antiparallel rotator, where we 
are able to confirm with our present studies many analytical predictions given by other 
authors in the past, in this paper for the first time the magnetospheres of the more 
general and complex system of inclinded and orthogonal rotators are investigated in 
the regime of the force-free surfaces (FFS) by the numerical calculation, using a full 
dynamical approach and taken into account radiation reaction. As in earlier work of 
our group, a 'standard set of parameters' is used. 

Under these conditions, the following results are found: Global charge separation 
exists for all degrees of inclination of the magnetic against the rotation axis with 
highest particle densities of lO^^cm^"'. Clouds of different charge are seperated by 
regions of vanishing particle number density. As expected, test particles inserted 
into the latter regions propagate into one of the adjacent clouds. The dependency 
of the electric monopole of the rotating sphere on the inclination angle is given by 
^« ~ IrT ^o^X- Furthermore strong polodial currents exist and locally averaged particle 
energies typically range up to 10"^^ — lO-*^^ eV, depending on the inclination angle. 

The results given by the presented work can be used as a starting point for an 
analytical description of neutron star magnetospheres in the near zone. A suitable 
approach have to consider a relativistic, non- neutral, charge seperated plasma. Although 
we prove corotation, an analytical approach should used a splitting in a corotational 
and non-corotational part of the description. In general, it is not usefull to assume a 
global force-free magnetosphere. The radiation reaction has to be taken into account, 
while using the ultra-relativistic approximation of the Lorentz-Dirac-Landau-equation 
is appropriate. 

Appendix A. Integration of the Continuity Equation 

This appendix describes the numerical integration of the continuity equation d^j^ = 
in detail. The used method is called flux corrected transport (FCT), and for details 
regarding general aspects of this method we refer to [38]. 

What follow we introduce a 'low-order' and a 'high-order' scheme in three spatial 
coordinates in order to construct a conservative FCT scheme. A spherical coordinate 
system is used and i,j,k denote the grid points to the e^, e^ and e(^-direction. 

The flux which flows from the cell i — 1 to cell i is called F, i u and the flux from 
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cell i to cell i+1 F-_^_ij^^. G and H are the corresponding fluxes in e^- and ei^-direction. 
With this notation a conservative discretisation is given by: 

^hj^k ^i,j,k 'Fi~{' 2 ,j ,k~^ '^i— 2 ,j ,k ^i,j-\-2,k~^^i,j—2,k '^i,j,k'j~2~^'^i,j,k~2 '("^'"^) 

As the 'low-order' scheme we use the Donor-Cell method in three dimensions 

= ^l,,k - 0.5(6 - 161) + 0.5(6 + 161) 

- 0.5(6 - 161) + 0.5(6 + 161) Ql,-,,k (A.2) 

- 0.5(6 - 161) QljMi + 0-5(6 + 161) (^7,,k-i ' 

with the Courant-numbers ^ concerning the three spatial coordinates which are indicated 
by the indices 1,2,3. The stability condition regarding this 'low-order' scheme is given 

by: 161 + 161 + 161 <i- 

A stable 'high-order' scheme (in three spatial coordinates) is developed by [40]: 

^n+i ^ ^"_At(V,(^v-0.5Atv(V,nv))). (A.3) 

Now wc arc able to write down the discretisation schemes. Using dimcnsionless 
units At' = LU At and Ar' = r^^ Ar and supressing the primes. The continuity equation 
is now given by dfN + diN — , where N is the particle density in the inertial frame 
of reference. 

The currents regarding the Donor-Cell method are given by 

{N /3^)j+i j^fe = - [ (/3r)i+ij,fe - \ {Pr)i+l,j,k\)^i+hj,k + 2 iiMi,j,k + \ {Mi,j,k\ Wij,k , 

{N(3^)- -^j = ^ [ (/3^)ij+i,fc - \ {(3^)ij+i^k\)Nij+i,k + ^ {{Pi))i,j,k + \ {Pi))i,j,k\ ]Nij,k , 

P'p)i,j,k+^ = 2 [ W'p)i,3,k+i ~ \ W'fi)i,j,k+l\)Nij^k+l + 2 iiP'p)i,j,k + \ Wip)i,j,k\ ]Nij,k ■ 

Using the following abbreviations for the surface elements Fl^, (i,j,k) = 2 sin -j^j sin A '^k , 
Fl^,(M,fc) = H^i+i sini^jAv'fc, F^(ij-fc) = -^^_i)^^^i the fluxes are 

given by: 

^i+i j,fc = Fl,. (,+ 1 [min( 0, (/3Ji+ij,fc ) iVi+ij,fc + max( 0, il3^)i,j,k ) iVi,i,fc] , 
Gij+i^k = Fltf, (ij+i,fe) [min( 0, {P^)i,j+i,k ) Nij+i^k + max( 0, {P^)i,j,k ) Nij^k] , 
Hij^k+^ = Fl^, (i,j,k+l) [min( 0, iP^)ij,j,k+i ) iV^j, fc+i + max( 0, {f3^)ij,k ) A^'ij.fc] . 

With these fluxes the 'low-order' scheme is given by 
1 

"^iyjjk 

Referring to (A.3) the 'high-order' scheme is given by 



^i',j',k - ^i,j,k + _ ^i+\,3,k + Fi-^,j,k Gi,j+^,k + Gi,j-\,k ^i,j,k+\ + ^i,j,k-i 



jyn+1 ^ - Ai 



ia. (r^ -^At ^)) + (sin^ -^At g)) 

rsmi? V 2 ^ y 
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with g := d,{N(3') = ±dr{r^N (3,) + -i^9^(sm^7V/3j + ^d^{N P^) . The fluxes 
regarding the 'high-order' scheme result in: 



+ 



+ 



— [sin^?,-+i( {Np^)ij+i,k + (A^/3^)i+ij+i,fe) 

4rj_,_i sm 

-sin^9,_i( (iV/5J„_i,fc + (iV/5^),+i,,_i,fc)] 

1 



4 1 sin 1? A 



r sin -i?^-^ 1 A 



(sin^ iV/3J. .+ i^,--At (sin^ 
[(sin^ A^/?^)ij+i,ik - (sin^ A^/3^)ij,d 



+ 



+ 



1 



4 Ar 



r^+i((7V/3,),+i,,-fc + (7V/3,),+ij+i,fc) 



-r^_,((iV/5,),_i,,-fc+(iV/3,), 



-ij+i 



1 



4rsini?. , 1 A 



ij+l,fe+l 

-((iV/?^)ij,fc-i + (iV/?^)ij+i,ik-i)] 



1 



+ 



+ 



r sin i9 A 
1 



4 Ar 
1 



4rsinT? A^? 



,k+l) 

-r^_i((iV/3,),_i,,,,+ (iV/5,),_i,,,fe+i); 



sin^,_i( (A^/3^),j_i,fc + 



)] 



Using the mentioned fluxes for the 'low-' and the 'high-order' scheme the 
construction of a full FCT scheme is straightforward (for details we refer to [38]). This 

2 2 2 

FCT scheme is stable if the following conditions are fulfilled: + 1^21^ + iCsl^ ^ 



l,At< 



V Ar y ~ Xrh.'ej ~\rsm^A'pJ 
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